A Counterfactual Simulation Model of Causal Judgment

1 Setup

1.1 Load packages

library("knitr")      # for knitting stuff
library("kableExtra") # for markdown tables
library("lme4")       # for linear mixed effects models
library("broom")      # for tidying model results 
library("brms")       # Bayesian regression with Stan
library("corrr")      # for tidy correlation matrix
library("xtable")     # for latex tables
library("jsonlite")   # for reading json files
library("janitor")    # cleans stuff
library("Hmisc")      # bootstrapped confidence intervals
library("tidybayes")  # for Bayesian data analysis
library("png")        # adding pngs to images
library("grid")       # functions for dealing with images 
library("egg")        # for geom_custom
library("patchwork")  # for figure panels
library("tidyverse")  # for wrangling, plotting, etc. 

1.2 Helper functions

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits) %>%
      print(include.rownames = F,
            booktabs = T)
  }
}

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>%
                drop_na(label),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme(text = element_text(size = 20),
          legend.position = "none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
          panel.spacing.y = unit(0.8, "cm"),
          plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
          axis.title.x = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% 
                 mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3,
                                               jitter.height = 0,
                                               dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme(text = element_text(size = 16),
          legend.position = "none",
          axis.title.x = element_blank(),
          panel.spacing.y = unit(c(.3), "cm"),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

# set ggplot theme 
theme_set(theme_classic())

2 Experiment 1: 2 balls

2.1 Clips

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(clip = 1:18,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 
                                             1, 1, 0, 1, 1, 1, 1, 1, 1),
                          outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 
                                                     0, 1, 1, 0, 0, 1, 0, 1, 1),
                          index_actual = rep(c("actual miss", 
                                               "actual close", 
                                               "actual hit"), each = 6),
                          index_counterfactual = rep(rep(c("counterfactual miss",
                                                           "counterfactual close", 
                                                           "counterfactual hit"),
                                                         each = 2),
                                                     3)) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp1.causal.long = df.exp1.causal$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

2.2.1 Demographics

# counterfactual condition 
df.exp1.counterfactual %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()

n	age_mean	age_sd	n_female	time_mean	time_sd
41	32	8.8	19	7.5	3.12

# causal condition 
df.exp1.causal %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()

n	age_mean	age_sd	n_female	time_mean	time_sd
41	33	10.3	21	9.65	6.58

2.3 Read in model predictions

2.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]

df.exp1.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp1.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
  group_by(clip) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>%
  ungroup() %>%
  left_join(df.exp1.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp1.counterfactual.model = df.exp1.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - 
                            df.exp1.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

2.3.2 Causal condition

df.exp1.causal.model = df.exp1.counterfactual.long %>% 
  group_by(clip) %>% 
  summarize(empirical = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.counterfactual.model %>% 
              select(-c(noise, sse)),
            by = "clip") %>% 
  rename(simulation = prediction) %>% 
  mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
         simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))

2.4 Counterfactual condition

2.4.1 Plots

set.seed(1)
df.plot = df.exp1.counterfactual.long %>% 
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.)/2)) %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", "actual close", "actual hit")),
         index_counterfactual = factor(index_counterfactual,
                                       levels = c("counterfactual miss",
                                                  "counterfactual close",
                                                  "counterfactual hit")),
         model = factor(model, levels = c("rating", "prediction"))) %>% 
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12


df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

2.4.2 Stats

2.4.2.1 Model evaluation

df.exp1.counterfactual.means %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise), 
            simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()

noise	simulation_r	simulation_rmse
1	0.97	10.95

2.5 Causal condition

2.5.1 Plots

set.seed(1)

model.name = c("empirical")
# model.name = c("simulation")

df.plot = df.exp1.causal.long %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp1.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, simulation, empirical),
               names_to = "model",
               values_to = "value") %>% 
  filter(model %in% c(model.name, "rating")) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index.actual = factor(index_actual, 
                          levels = c("actual miss", "actual close", "actual hit")),
    index.counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit"))) %>%
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)

# ggsave("../../figures/plots/exp1_causal_bars.pdf",
#        width = 8,
#        height = 6)

2.5.2 Stats

2.5.2.1 Model evaluation

df.exp1.causal.long %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, empirical),
            empirical_rmse = rmse(rating, empirical),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  print_table()

empirical_r	empirical_rmse	simulation_r	simulation_rmse
0.96	8.57	0.93	15.15

2.5.2.2 Bayesian regression

# random intercepts model 
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + 
                            outcome_actual + (1 + index_actual + index_counterfactual + 
                                                outcome_actual | participant),
                          data = df.exp1.causal.long,
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp1_causal")

2.5.2.3 Specific hypotheses

df.exp1.causal.posterior = df.exp1.causal.long %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp1_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

func_posterior_difference = function(data, clip1, clip2){
  data %>% 
    mutate(difference = .data[[clip1]] - .data[[clip2]]) %>% 
    pull(difference) %>% 
    mean_hdci() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>% 
    print()  
}

func_posterior_difference(data = df.exp1.causal.posterior,
                          clip1 = "8",
                          clip2 = "13")

            difference
1 (0.26 [-6.51, 7.31])

2.5.3 Tables

fit.brm_exp1_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()

term	estimate	std.error	lower 95% HDI	upper 95% HDI
intercept	-18.33	3.36	-23.81	-12.78
index_actualactualclose	4.31	3.40	-1.40	9.83
index_actualactualhit	4.06	4.84	-4.01	11.88
index_counterfactualcounterfactualclose	-17.52	3.37	-22.96	-12.01
index_counterfactualcounterfactualhit	-61.37	4.38	-68.54	-54.30
outcome_actual	92.36	5.11	84.04	100.96

3 Experiment 2: Brick & Teleport

3.1 Clips

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
                                header = T,
                                stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
                                        header = T,
                                        stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(clip = 1:20,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
                                             0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
                          outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 
                                                     1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
                          index_actual = c(rep(c("actual miss", 
                                                 "actual close",
                                                 "actual hit"),
                                               each = 6),
                                           rep("practice", 2)),
                          index_counterfactual = c(rep(rep(c("counterfactual miss",
                                                             "counterfactual close",
                                                             "counterfactual hit"),
                                                           each = 2), 
                                                       3),
                                                   rep("practice", 2))) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp2.causal_first.long = df.exp2.causal_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
         condition = "causal_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
         condition = "counterfactual_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + 
                       max(df.exp2.causal_first.long$participant)))

3.2.1 Demographics

# counterfactual condition 
df.exp2.counterfactual_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time) %>% 
  bind_rows(df.exp2.causal_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time)) %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()

n	age_mean	age_sd	n_female	time_mean	time_sd
82	35	12.2	45	21.25	5.11

3.3 Read in model predictions

3.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]

df.exp2.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp2.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp2.counterfactual.means = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  # mutate(rating = abs(rating)) %>% 
  group_by(clip) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp2.counterfactual.model = df.exp2.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

3.3.2 Causal condition

# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
  filter(question == "counterfactual",
    clip <= 18) %>%
  group_by(clip, condition) %>%
  summarize(rating = mean(rating)) %>%
  pivot_wider(names_from = condition,
              values_from = rating) %>% 
  left_join(df.exp2.combined.long %>%
              filter(question == "counterfactual", clip <= 18) %>%
              group_by(clip) %>%
              summarize(combined = mean(rating)),
            by = "clip") %>% 
  left_join(df.exp2.counterfactual.model,
            by = "clip") %>% 
  select(-c(sse, noise)) %>% 
  rename(simulation = prediction) %>% 
  mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
            .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))

3.4 Counterfactual condition

3.4.1 Plots

set.seed(1)

df.plot = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.counterfactual.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")),
         model = factor(model,
                             levels = c("rating", "prediction")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

3.4.2 Stats

3.4.2.1 Model evaluation

df.exp2.counterfactual.means %>% 
  left_join(df.exp2.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise),
            simulation_r = cor(rating, prediction),
            simulation_rmse = rmse(rating, prediction)) %>% 
  print_table()

noise	simulation_r	simulation_rmse
0.9	0.96	13.46

3.4.2.2 Bayesian regressions

3.4.2.2.1 Check for potential order effects

# random intercepts & slopes model 
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * 
                                    (index_actual + index_counterfactual) + 
                                    (1 + index_actual + 
                                       index_counterfactual | participant),
                                  data = df.exp2.combined.long %>% 
                                    na.omit() %>% 
                                    filter(question == "counterfactual"),
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp2_counterfactual")

3.4.3 Tables

fit.brm_exp2_counterfactual %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  print_table()

term	estimate	std.error	lower 95% HDI	upper 95% HDI
intercept	13.38	2.58	9.12	17.55
conditioncounterfactual_first	-1.11	3.65	-6.99	4.82
index_actualactualclose	-1.48	2.78	-6.02	3.11
index_actualactualhit	-4.42	2.40	-8.32	-0.46
index_counterfactualcounterfactualclose	37.94	3.33	32.39	43.35
index_counterfactualcounterfactualhit	73.80	4.17	67.00	80.57
conditioncounterfactual_first:index_actualactualclose	-0.55	3.84	-7.11	5.73
conditioncounterfactual_first:index_actualactualhit	-0.85	3.41	-6.52	4.85
conditioncounterfactual_first:index_counterfactualcounterfactualclose	1.23	4.79	-6.62	9.18
conditioncounterfactual_first:index_counterfactualcounterfactualhit	-0.02	5.92	-9.86	9.55

  # xtable() %>%
  # print(include.rownames = F,
  #       booktabs = T)

3.5 Causal condition

3.5.1 Plots

set.seed(1)

func_exp2_causal_plot = function(condition.name, model.name){

df.plot = df.exp2.combined.long %>%
  filter(question == "causal", clip <= 18) %>%
  filter(condition %in% condition.name) %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c("rating", all_of(model.name)),
               names_to = "model",
               values_to = "value") %>% 
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual, 
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")
}

plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
                                       model.name = "combined"),
                 func_exp2_causal_plot(condition.name = "causal_first",
                                       model.name = "combined"))

# model.name = "combined"
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# 
# func_exp2_causal_plot(condition.name = condition.name,
#                       model.name = model.name),
# 
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)

3.5.2 Stats

3.5.2.1 Model evaluation

df.data = df.exp2.combined.long %>%
  filter(question == "causal",
         clip <= 18) %>%
  mutate(rating = abs(rating)) %>%
  group_by(condition,
           clip,
           outcome_actual,
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual"))

df.data %>% 
  group_by(condition) %>% 
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined)) %>% 
  print_table()

condition	empirical_r	empirical_rmse
causal_first	0.92	13.78
counterfactual_first	0.99	5.27

df.data %>%
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>%
  print_table()

empirical_r	empirical_rmse	simulation_r	simulation_rmse
0.95	10.43	0.92	19.09

3.5.2.2 Bayesian regression

# random intercepts + slopes model 
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * 
                            (index_actual + index_counterfactual + outcome_actual) +
                            (1 + index_actual + index_counterfactual + 
                               outcome_actual | participant),
                          data = df.exp2.combined.long %>% 
                            filter(question == "causal",
                                   clip <= 18),
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp2_causal")

3.5.2.3 Specific hypotheses

df.exp2.causal.posterior = df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp2_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, condition, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

# difference between clips 13 and 14 versus 8 
df.exp2.causal.posterior %>% 
    mutate(difference = (`13`  + `14`) / 2 - `8`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

               difference
1 (-0.55 [-12.09, 10.97])

# difference between clips 5 and 6 versus 11 
df.exp2.causal.posterior %>% 
    mutate(difference = (`5`  + `6`) / 2 - `11`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

             difference
1 (-5.48 [-13.8, 2.73])

3.5.3 Tables

3.5.3.1 Bayesian regression

fit.brm_exp2_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()

term	estimate	std.error	lower 95% HDI	upper 95% HDI
intercept	-24.07	3.50	-29.80	-18.27
conditioncounterfactual_first	11.47	4.91	3.29	19.40
index_actualactualclose	8.17	3.44	2.39	13.62
index_actualactualhit	13.42	4.90	5.42	21.30
index_counterfactualcounterfactualclose	-30.20	3.72	-36.41	-24.36
index_counterfactualcounterfactualhit	-50.29	4.60	-57.94	-42.84
outcome_actual	98.53	5.40	89.53	107.57
conditioncounterfactual_first:index_actualactualclose	-5.39	4.94	-13.38	2.80
conditioncounterfactual_first:index_actualactualhit	-16.97	7.02	-28.61	-5.75
conditioncounterfactual_first:index_counterfactualcounterfactualclose	-8.13	5.17	-16.51	0.35
conditioncounterfactual_first:index_counterfactualcounterfactualhit	-25.51	6.44	-36.53	-15.18
conditioncounterfactual_first:outcome_actual	10.72	7.60	-1.63	23.34

4 Experiment 3: 3 balls

4.1 Clips

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
                          outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
                                           0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
                          outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                          outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
                          outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                          clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", "condition"),
                    n()/6),
         participant = rep(1:(n()/6), each = 6)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "rating"), n()/2),
         participant = rep(1:(n()/64), each = 64),
         order = rep(rep(1:32, each = 2), n()/64)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  left_join(df.exp3.counterfactual.demographics %>% 
              select(participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>% 
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
                          header = T,
                          stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", 
                      "counterbalance", "replayType", "experiment"),
                    n()/8),
         participant = rep(1:(n()/8), each = 8)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)

# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
         participant = rep(1:(n()/(32*6)), each = (32*6)),
         order = rep(rep(1:32, each = 6), n()/(32*6))) %>% 
  filter(!name %in% c("x", "y", "z")) %>% 
pivot_wider(names_from = name,
            values_from = value) %>% 
  left_join(df.exp3.causal.demographics %>% 
              select(participant, counterbalance),
            by = "participant") %>% 
  mutate(A = ifelse(counterbalance == 1, rating1, rating2),
         B = ifelse(counterbalance == 1, rating2, rating1)) %>% 
  select(-rating1, -rating2) %>% 
  pivot_longer(cols = c(A, B),
               names_to = "ball",
               values_to = "rating") %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.2.1 Demographics

# counterfactual condition 
df.exp3.counterfactual.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()

n	age_mean	age_sd	n_female	time_mean	time_sd
80	33	10.1	34	18.08	4.63

# causal condition 
df.exp3.causal.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()

n	age_mean	age_sd	n_female	time_mean	time_sd
41	34	10.5	21	21.19	4.96

4.3 Read in model predictions

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]

df.exp3.model = files %>%
  map_dfr(~ read_csv(str_c(path, .))) %>% 
  rename(clip = trial)

4.3.1 Counterfactual condition

# calculate mean counterfactual judgments 
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp3.model.counterfactual = df.exp3.model %>% 
  group_by(noise) %>% 
  select(clip, contains("whether"), noise) %>% 
  pivot_longer(cols = c(A_whether, B_whether),
               names_to = "ball",
               values_to = "prediction") %>% 
  mutate(ball = str_remove(ball, "_whether")) %>% 
  arrange(noise, clip, ball) %>% 
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data,
                   ~ sum((.x$prediction*100 -
                            df.exp3.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

4.3.2 Causal condition

# calculate mean causal judgments 
df.exp3.causal.means = df.exp3.causal.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

df.exp3.model.causal = df.exp3.model %>% 
  # take into account difference-making 
  mutate_at(.vars = vars(A_how:A_robust),
            .funs = list(~ . * A_difference)) %>% 
  mutate_at(.vars = vars(B_how:B_robust),
            .funs = list(~ . * B_difference)) %>% 
  # choose model based on best fit with counterfactual data 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>% 
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(clip, ball:robust)

4.3.3 Heuristic model

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
             "ballA_velx",
             "ballA_vely",
             "ballB_velx",
             "ballB_vely",
             "ballE_velx",
             "ballE_vely")) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
            mutate(clip = i) %>%
            select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(ball == "ballE",
         str_detect(object, "goal|Left|Right")) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp, by = "clip") %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  count(clip, ball, object) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
  group_by(clip) %>%
  summarize(AE = any(AE %>% as.logical()) * 1,
            BE = any(BE %>% as.logical()) * 1) %>%
  mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
         B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
  select(clip, A, B) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "exclusive")

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ball"),
         ball == "ballE") %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(E.speed.diff = speed.diff,
         E.direction.diff = direction.diff)

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(total.speed.diff = speed.diff,
         total.direction.diff = direction.diff)

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(str_detect(ball, "ballA|ballB"),
         str_detect(object, "ball")) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
         AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarize(A = any(!is.na(AE)),
            B = any(!is.na(BE)),
            A = ifelse(any(!is.na(AB)) &
                         any(!is.na(BE)) &
                         max(AB, na.rm = T) < min(BE, na.rm = T),
                       T,
                       A),
            B = ifelse(any(!is.na(AB)) &
                         any(!is.na(AE)) &
                         max(AB, na.rm = T) < min(AE, na.rm = T),
                       T,
                       B)) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "transfer") %>% 
  arrange(clip, ball)

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
         speed = abs(velx) + abs(vely)) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.4.1 Plots

4.4.1.1 Bar graph

set.seed(1)

df.plot = df.exp3.counterfactual.long %>%
  left_join(df.exp3.model.counterfactual %>%
              select(clip, ball, prediction) %>%
              mutate(ball = as.factor(ball)),
            by = c("clip", "ball")) %>%
  left_join(df.exp3.clipinfo,
            by = "clip") %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "prediction", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "prediction"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>%
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
#        width = 12,
#        height = 6)

4.4.2 Stats

4.4.2.1 Model evaluation

# best fitting model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.model.counterfactual,
            by = c("clip", "ball")) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()

simulation_r	simulation_rmse
0.86	19.84

# deterministic model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_a, outcome_b) %>% 
              pivot_longer(cols = -clip,
                           names_to = "ball",
                           values_to = "outcome") %>% 
              mutate(ball = factor(ball,
                                   levels = c("outcome_a", "outcome_b"),
                                   labels = c("B", "A"))),
            by = c("clip", "ball")) %>% 
  mutate(outcome = outcome * 100) %>% 
  summarize(simulation_r = cor(rating_mean, outcome),
            simulation_rmse = rmse(rating_mean, outcome)) %>% 
  print_table()

simulation_r	simulation_rmse
0.83	30.28

4.5 Causal condition

4.5.1 Stats

4.5.1.1 Bayesian mixed effects model

df.data = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# using whether-causation from the model 
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
                                  data = df.data,
                                  save_all_pars = T,
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp3_causal_w")

fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
                               (1 + whether + how | participant),
                                   data = df.data,
                                   save_all_pars = T,
                                   cores = 2,
                                   seed = 1,
                                   file = "cache/brm_exp3_causal_wh")

fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
                                (1 + whether + how + sufficient | participant),
                                    data = df.data,
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_whs")

# fit.brm_exp3_causal_w %>% summary()
# fit.brm_exp3_causal_wh %>% summary()
# fit.brm_exp3_causal_whs %>% summary()

4.5.1.2 Model comparison

fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
                                            criterion = c("loo", "waic"),
                                            reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
                                             criterion = c("loo", "waic"),
                                             reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
                                              criterion = c("loo", "waic"),
                                              reloo = T)

No problematic observations found. Returning the original 'loo' object.

loo_compare(fit.brm_exp3_causal_w,
            fit.brm_exp3_causal_wh,
            fit.brm_exp3_causal_whs)

                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

• z-scored predictors and invidiual responses

# Regression based on features
df.regression.features = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  left_join(df.features %>% 
              mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
            by = c("clip", "ball")) %>%
  clean_names() %>% 
  mutate(e_moving = 1 - e_moving)

# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving + 
                                      speed + 
                                      contact_e + 
                                      e_speed_diff + 
                                      e_direction_diff + 
                                      total_speed_diff + 
                                      total_direction_diff + 
                                      transfer + 
                                      e_moving + 
                                      exclusive + (1 | participant),
                                    prior = priors,
                                    data = df.regression.features %>% 
                                      select(-c(clip, ball, order, clipindex,
                                                contains("outcome"))),
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.4 Predictions for mean causal judgments

func_fit_data = function(fit){
  fit %>% 
    fitted(newdata = df.exp3.model.causal %>% 
             left_join(df.features %>% 
                         mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
                       by = c("clip", "ball")) %>%
             clean_names() %>% 
             mutate(e_moving = 1 - e_moving),
           re_formula = NA) %>% 
    as_tibble() %>% 
    pull(Estimate)
}

df.exp3.causal.regression = df.exp3.causal.means %>% 
  mutate(w = func_fit_data(fit.brm_exp3_causal_w),
         wh = func_fit_data(fit.brm_exp3_causal_wh),
         whs = func_fit_data(fit.brm_exp3_causal_whs),
         heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))

4.5.1.5 Model evaluation

prediction = "whs"

df.exp3.causal.regression %>% 
  summarize(simulation_r = cor(rating_mean, .[[prediction]]),
            simulation_rmse = rmse(rating_mean, .[[prediction]])) %>% 
  print_table()

simulation_r	simulation_rmse
0.87	13.07

4.5.1.6 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = 
  brm(formula = rating ~ 1,
      data = df.fit %>% 
        filter(participant == 1),
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = 
  brm(formula = rating ~ 1 + whether,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = 
  brm(formula = rating ~ 1 + whether + how,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = 
  brm(formula = rating ~ 1 + whether + how + sufficient,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ rmse(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ rmse(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ rmse(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "data/exp3_causal_individual_fit.RData")

4.5.1.7 Load individual participant regressions

load("data/exp3_causal_individual_fit.RData")

# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>% 
  count(best_model) %>% 
  print_table()

best_model	n
wh	2
whs	39

4.5.2 Plots

4.5.2.1 Bar graph

set.seed(1)
model_index = "whs"

# model predictions 
model_prediction = list(fit.brm_exp3_causal_w,
                        fit.brm_exp3_causal_wh,
                        fit.brm_exp3_causal_whs) %>%
  map_dfr(~ fitted(., df.exp3.causal.means %>% 
                     left_join(df.exp3.model.causal,
                               by = c("clip", "ball")),
                   re_formula = NA) %>% 
            as_tibble()) %>% 
  mutate(ball = rep(c("A", "B"), n()/2),
         clip = rep(rep(1:32, each = 2), 3),
         model = rep(c("w", "wh", "whs"),
                     each = 64))

df.plot = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_both),
            by = c("clip")) %>% 
  left_join(model_prediction %>% 
              filter(model == model_index) %>% 
              select(model = Estimate, clip, ball),
            by = c("clip", "ball")) %>% 
  pivot_longer(cols = c(rating, model),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
         colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
         colorindex = ifelse(model == "model", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "model"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>% 
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)

# ggsave("../../figures/plots/exp3_causal_bars.pdf",
#        width = 12,
#        height = 6)

4.5.2.2 Selection

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarize(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  pivot_longer(cols = c(mean, w, wh, whs),
               names_to = "index",
               values_to = "value") %>% 
  mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain", "double prevention", "joint causation",
                                  "overdetermination", "preemption")))

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

func_scatterplot = function(model){
  if(model == "heuristic"){
    xlabel = "Heuristic"
  }else{
    xlabel = bquote(CSM[.(toupper(model))])
  }
  
  l.models = list(w = fit.brm_exp3_causal_w, 
                  wh = fit.brm_exp3_causal_wh,
                  whs = fit.brm_exp3_causal_whs,
                  heuristic = fit.brm_exp3_causal_heuristic)
  
  df.data = df.exp3.causal.means %>% 
    left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
    left_join(df.regression.features, by = c("clip", "ball"))
  
  df.plot = l.models[[model]] %>%
    fitted(newdata = df.data,
           re_formula = NA) %>% 
    as_tibble() %>% 
    clean_names() %>%
    bind_cols(df.data)
  
  p = ggplot(data = df.plot,
             mapping = aes(x = estimate, 
                           y = rating_mean)) +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) + 
    geom_smooth(aes(y = estimate,
                    ymin = q2_5,
                    ymax = q97_5),
                stat = "identity",
                color = "black") +
    geom_linerange(size = 0.75, 
                   mapping = aes(ymin = rating_low,
                                 ymax = rating_high),
                   color = "gray50") +
    geom_point(size = 2) + 
    scale_color_manual(values = c("black", "#e41a1c"),
                       guide = F) +
    coord_fixed(xlim = c(0, 100),
                ylim = c(0, 100)) +
    labs(x = xlabel,
         y = "responsibility judgment") +
    annotate(geom = "text",
             label = str_c(
               "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2) %>% 
                 as.character() %>% 
                 str_sub(start = 2),
               "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2)),
             x = 5,
             y = 95,
             size = 6,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 100, 25)) +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    theme(text = element_text(size = 20))
  return(p)
}

# for creating and saving an individual scatter plot 
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"

# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

list(func_scatterplot(model = "w"),
     func_scatterplot(model = "wh"),
     func_scatterplot(model = "whs"),
     func_scatterplot(model = "heuristic")) %>%
  wrap_plots(ncol = 2)

4.5.2.4 Individual participant variance for a selection of clips (clustered)

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.cluster = df.exp3.causal.individual_fit %>%
  select(participant, fit_whs) %>%
  mutate(estimates = map(fit_whs, tidy)) %>%
  select(participant, estimates) %>%
  unnest(estimates) %>%
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>%
  mutate(term = str_remove(term, "b_")) %>%
  select(participant, term, estimate) %>%
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate(group = ifelse(how > whether, "how", "whether"))

df.plot = df.plot %>% 
  left_join(df.cluster %>% 
              select(participant, group),
            by = "participant")

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

# ggplot(df.plot,aes(x=ball,y=rating,group=id,color = best))+
p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(group == "how"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(group == "whether"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#e41a1c", "#377eb8"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

print(p)

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.2.5 Ternary plots for individual regression results

• doesn’t work at the moment: ggtern is incompatible with the most recent ggplot2 update

library("ggtern")
df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, fit_whs) %>% 
  mutate(estimates = map(fit_whs, tidy)) %>% 
  select(participant, estimates) %>% 
  unnest(estimates) %>% 
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>% 
  mutate(term = str_remove(term, "b_")) %>% 
  select(participant, term, estimate) %>% 
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate_at(vars(how:whether),
            .funs = list(norm = ~ . / (how + whether + sufficient))) %>% 
  mutate(color = 0,
         color = ifelse(how_norm == max(how_norm), 1, color),
         color = ifelse(whether_norm == max(whether_norm), 2, color),
         color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
         color = factor(color))

ggplot(data = df.plot,
       mapping = aes(x = how,
                     y = sufficient,
                     z = whether,
                     color = color)) + 
  geom_point(alpha = 0.7,
             size = 2,
             show.legend = F) + 
  scale_color_manual(values = c("black", "red", "green", "blue")) + 
  coord_tern() +
  # theme_bw() + 
  theme_showarrows() + 
  theme(text = element_text(size = 20),
                tern.axis.title.T = element_blank(),
                tern.axis.title.L = element_blank(),
                tern.axis.title.R = element_blank())

# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
#        width = 5,
#        height = 5)

4.5.3 Tables

4.5.3.1 Relationship between predictors

df.exp3.model %>% 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(difference, whether, how, sufficient, robust) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  # print_table(format = "latex")
  print_table()

rowname	difference	whether	how	sufficient	robust
difference
whether	.50
how	.79	.27
sufficient	.21	.10	.36
robust	.43	.93	.24	.20

4.5.3.2 Population level predictors in the mixed effects models

fit.brm_exp3_causal_w %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(model = "CSM_w") %>% 
  bind_rows(fit.brm_exp3_causal_wh %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_wh")) %>% 
  bind_rows(fit.brm_exp3_causal_whs %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_whs")) %>% 
  mutate(term = tolower(term),
         term = str_remove_all(term, "b_")) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  select(model, everything()) %>% 
  # print_table(format = "latex")
  print_table()

model	term	estimate	std.error	lower 95% HDI	upper 95% HDI
CSM_w	intercept	31.68	1.61	29.02	34.37
CSM_w	whether	39.30	2.16	35.76	42.82
CSM_wh	intercept	10.29	1.56	7.80	12.83
CSM_wh	whether	28.40	2.68	23.96	32.70
CSM_wh	how	36.07	2.60	31.80	40.27
CSM_whs	intercept	10.43	1.59	7.82	13.01
CSM_whs	whether	27.93	2.60	23.67	32.18
CSM_whs	how	26.47	2.91	21.76	31.35
CSM_whs	sufficient	30.90	3.06	25.77	35.96

4.5.3.3 Individual participants regression results

df.exp3.causal.individual_fit %>% 
  select_if(~ (is.numeric(.) | is.factor(.))) %>% 
  mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>% 
  select(participant, everything(), best_model) %>% 
  # print_table(format = "latex")
  print_table()

participant	r_w	r_wh	r_whs	rmse_w	rmse_wh	rmse_whs	best_model
1	0.29	0.38	0.48	30.00	29.03	27.74	whs
2	0.20	0.77	0.77	39.32	27.94	27.55	whs
3	0.37	0.78	0.79	38.90	28.00	26.95	whs
4	0.41	0.55	0.63	43.62	40.60	38.31	whs
5	0.45	0.53	0.54	39.23	37.32	36.81	whs
6	0.18	0.54	0.54	40.84	36.32	36.07	whs
7	0.49	0.61	0.64	32.96	29.93	29.15	whs
8	0.39	0.63	0.68	38.37	33.26	31.64	whs
9	0.40	0.72	0.76	35.86	28.34	26.39	whs
10	0.28	0.52	0.53	29.13	26.28	25.85	whs
11	0.51	0.56	0.60	33.66	32.19	31.17	whs
12	0.56	0.57	0.70	32.78	32.08	28.84	whs
13	0.56	0.71	0.74	29.83	25.40	24.21	whs
14	0.43	0.47	0.50	33.72	32.82	32.23	whs
15	0.19	0.35	0.37	39.77	38.29	37.94	whs
16	0.62	0.67	0.76	40.49	37.98	34.39	whs
17	0.42	0.66	0.71	28.57	23.82	22.42	whs
18	0.31	0.35	0.39	37.78	37.09	36.57	whs
19	0.37	0.60	0.62	33.41	29.23	28.62	whs
20	0.57	0.61	0.68	36.47	35.00	32.71	whs
21	0.43	0.50	0.57	30.95	29.61	28.32	whs
22	0.51	0.67	0.71	33.72	29.36	28.00	whs
23	0.33	0.45	0.52	38.76	37.00	35.57	whs
24	0.46	0.63	0.69	34.01	29.96	28.22	whs
25	0.61	0.66	0.70	30.60	28.91	27.55	whs
26	0.57	0.71	0.76	40.43	35.03	32.80	whs
27	0.46	0.52	0.66	43.12	41.49	38.12	whs
28	0.30	0.45	0.46	29.33	27.62	27.38	whs
29	0.48	0.60	0.65	38.93	36.03	34.48	whs
30	0.18	0.70	0.70	33.06	24.95	24.71	whs
31	0.48	0.59	0.69	37.54	34.72	31.76	whs
32	0.37	0.61	0.70	40.62	35.94	32.92	whs
33	0.50	0.56	0.64	38.64	36.79	34.85	whs
34	0.24	0.51	0.52	44.50	40.86	40.19	whs
35	0.15	0.32	0.30	34.56	33.36	33.37	wh
36	0.19	0.84	0.82	45.38	29.58	29.70	wh
37	0.55	0.63	0.67	36.04	33.36	32.11	whs
38	0.52	0.61	0.67	36.03	33.39	31.53	whs
39	0.47	0.77	0.77	29.61	21.75	21.58	whs
40	0.46	0.78	0.79	32.41	23.60	22.69	whs
41	0.32	0.50	0.58	32.90	30.38	28.69	whs

4.5.3.4 CSMwhs predictions for selection of cases

df.exp3.model.causal %>%
  left_join(df.exp3.causal.regression,
            by = c("clip", "ball")) %>% 
  filter(clip %in% c(7, 23, 16, 3, 15)) %>%
  mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
            .funs = ~ round(., 2)) %>%
  select(clip, ball, difference, whether, how, sufficient, robust) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ ifelse(. < 0.5,
                             str_c("xmark (", ., ")"),
                             str_c("cmark (", ., ")"))) %>%
  arrange(clip) %>% 
  # print_table(format = "latex")
  print_table()

clip	ball	difference	whether	how	sufficient	robust
7	A	cmark (1)	xmark (0.34)	cmark (1)	xmark (0)	xmark (0.25)
7	B	cmark (1)	cmark (1)	cmark (1)	cmark (0.67)	cmark (0.6)
23	A	xmark (0.05)	xmark (0)	xmark (0)	xmark (0)	xmark (0)
23	B	cmark (0.91)	cmark (0.79)	xmark (0)	xmark (0)	cmark (0.72)
16	A	cmark (1)	xmark (0.23)	cmark (1)	cmark (1)	xmark (0.35)
16	B	xmark (0)	xmark (0)	xmark (0)	xmark (0)	xmark (0)
3	A	cmark (1)	cmark (0.88)	cmark (1)	xmark (0.12)	cmark (0.76)
3	B	cmark (1)	cmark (0.89)	cmark (1)	xmark (0.11)	cmark (0.75)
15	A	cmark (1)	xmark (0.01)	cmark (1)	cmark (0.99)	xmark (0.1)
15	B	cmark (1)	xmark (0.01)	cmark (1)	cmark (1)	xmark (0.1)

4.5.3.5 CSMwhs predictions for all cases

df.exp3.causal.regression %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball")) %>% 
  left_join(df.exp3.clipinfo %>% 
              select(-clipindex),
            by = c("clip")) %>% 
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ . * 100) %>% 
  mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
  select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust,
         w, wh, whs, heuristic, rating = rating_mean) %>% 
  # print_table(format = "latex", digits = 0)
  print_table(digits = 0)

clip	ball	outcome_both	outcome_a	outcome_b	outcome_none	difference	whether	how	sufficient	robust	w	wh	whs	heuristic	rating
1	A	0	0	0	0	100	40	100	23	36	47	58	55	57	42
1	B	0	0	0	0	100	15	100	16	9	38	51	46	54	37
2	A	0	0	0	0	57	12	0	0	10	36	14	14	25	21
2	B	0	0	0	0	18	0	0	0	0	32	10	10	24	19
3	A	1	0	0	0	100	88	100	12	76	66	71	65	72	76
3	B	1	0	0	0	100	89	100	11	75	67	72	65	72	75
4	A	1	0	0	0	100	78	100	4	78	62	69	60	58	63
4	B	1	0	0	0	100	95	100	15	57	69	73	68	54	78
5	A	0	0	1	0	100	90	100	0	47	67	72	62	47	71
5	B	0	0	1	0	100	0	100	0	0	32	46	37	68	22
6	A	0	0	1	0	100	59	100	16	35	55	63	59	53	73
6	B	0	0	1	0	100	18	100	6	14	39	52	44	53	22
7	A	1	0	1	0	100	34	100	0	25	45	56	46	70	59
7	B	1	0	1	0	100	100	100	67	60	71	75	86	64	79
8	A	1	0	1	0	0	0	0	0	0	32	10	10	25	7
8	B	1	0	1	0	100	100	100	100	100	71	75	96	84	92
9	A	0	1	0	0	0	0	0	0	0	32	10	10	14	8
9	B	0	1	0	0	100	100	100	0	100	71	75	65	78	90
10	A	0	1	0	0	77	18	0	0	22	39	15	16	15	23
10	B	0	1	0	0	98	79	0	0	63	63	33	33	21	55
11	A	1	1	0	0	100	70	100	77	68	59	66	80	71	93
11	B	1	1	0	0	0	0	0	0	0	32	10	10	16	4
12	A	1	1	0	0	100	82	100	74	83	64	70	83	57	77
12	B	1	1	0	0	100	0	100	12	24	32	46	41	53	37
13	A	0	1	1	0	67	34	0	0	35	45	20	20	14	8
13	B	0	1	1	0	70	35	0	0	35	46	20	20	21	64
14	A	0	1	1	0	97	91	0	0	59	67	36	36	21	22
14	B	0	1	1	0	91	77	0	0	51	62	32	32	20	18
15	A	1	1	1	0	100	1	100	99	10	32	47	68	71	76
15	B	1	1	1	0	100	1	100	100	10	32	47	68	71	76
16	A	1	1	1	0	100	23	100	100	35	41	53	74	80	92
16	B	1	1	1	0	0	0	0	0	0	32	10	10	14	4
17	A	0	0	0	1	100	19	100	37	18	39	52	54	66	69
17	B	0	0	0	1	100	0	100	36	17	32	46	48	65	46
18	A	0	0	0	1	100	11	100	40	17	36	49	52	55	63
18	B	0	0	0	1	100	7	100	37	9	34	48	50	56	66
19	A	1	0	0	1	100	74	100	7	65	61	67	60	55	53
19	B	1	0	0	1	100	72	100	7	65	60	67	59	55	49
20	A	1	0	0	1	100	92	100	8	72	68	72	65	57	41
20	B	1	0	0	1	100	88	100	4	53	66	71	63	56	71
21	A	0	0	1	1	100	47	100	40	45	50	60	63	58	80
21	B	0	0	1	1	100	9	100	21	10	35	49	46	59	18
22	A	0	0	1	1	100	100	100	89	83	71	75	92	48	60
22	B	0	0	1	1	100	8	100	0	15	35	49	39	53	42
23	A	1	0	1	1	5	0	0	0	0	32	10	10	15	3
23	B	1	0	1	1	91	79	0	0	72	63	33	33	22	39
24	A	1	0	1	1	100	66	100	4	63	58	65	57	57	44
24	B	1	0	1	1	100	94	100	22	79	69	73	70	54	73
25	A	0	1	0	1	100	25	100	21	26	42	53	50	69	43
25	B	0	1	0	1	100	74	100	54	65	61	67	74	56	73
26	A	0	1	0	1	100	6	100	3	9	34	48	40	60	39
26	B	0	1	0	1	100	87	100	35	54	66	71	72	46	69
27	A	1	1	0	1	100	97	100	52	97	70	74	80	67	80
27	B	1	1	0	1	0	0	0	0	0	32	10	10	17	6
28	A	1	1	0	1	100	90	100	22	80	67	72	69	79	89
28	B	1	1	0	1	0	0	0	0	0	32	10	10	12	5
29	A	0	1	1	1	100	58	100	24	44	55	63	61	66	47
29	B	0	1	1	1	100	63	100	24	38	57	64	62	54	67
30	A	0	1	1	1	100	57	100	29	49	54	63	62	63	58
30	B	0	1	1	1	100	46	100	24	39	50	60	57	63	56
31	A	1	1	1	1	100	2	100	4	3	32	47	39	63	44
31	B	1	1	1	1	100	4	100	4	4	33	47	39	53	46
32	A	1	1	1	1	0	0	0	0	0	32	10	10	16	5
32	B	1	1	1	1	100	75	100	66	73	61	68	78	65	71

4.5.3.6 Heuristic model

fit.brm_exp3_causal_heuristic %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()

term	estimate	std.error	lower 95% HDI	upper 95% HDI
intercept	49.79	1.46	47.39	52.23
moving	0.22	0.22	0.01	0.66
speed	2.06	0.85	0.64	3.43
contact_e	0.38	0.35	0.02	1.07
e_speed_diff	0.12	0.12	0.01	0.37
e_direction_diff	1.04	0.74	0.08	2.43
total_speed_diff	2.21	0.95	0.67	3.80
total_direction_diff	3.99	0.91	2.45	5.45
transfer	15.59	0.79	14.30	16.90
e_moving	0.18	0.17	0.01	0.51
exclusive	4.39	0.73	3.16	5.56

5 Session info

sessionInfo()

R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.0    stringr_1.4.0    dplyr_0.8.5      purrr_0.3.3     
 [5] readr_1.3.1      tidyr_1.0.2      tibble_2.1.3     tidyverse_1.3.0 
 [9] patchwork_1.0.0  egg_0.4.5        gridExtra_2.3    png_0.1-7       
[13] tidybayes_2.0.2  Hmisc_4.3-1      ggplot2_3.3.0    Formula_1.2-3   
[17] survival_3.1-11  lattice_0.20-40  janitor_1.2.1    jsonlite_1.6.1  
[21] xtable_1.8-4     corrr_0.4.2      brms_2.12.0      Rcpp_1.0.4      
[25] broom_0.5.5      lme4_1.1-21      Matrix_1.2-18    kableExtra_1.1.0
[29] knitr_1.28      

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.1.5      plyr_1.8.6          
  [4] igraph_1.2.5         splines_3.6.2        svUnit_0.7-12       
  [7] crosstalk_1.1.0.1    rstantools_2.0.0     inline_0.3.15       
 [10] digest_0.6.25        htmltools_0.4.0      rsconnect_0.8.16    
 [13] fansi_0.4.1          magrittr_1.5         checkmate_2.0.0     
 [16] cluster_2.1.0        modelr_0.1.6         matrixStats_0.56.0  
 [19] xts_0.12-0           prettyunits_1.1.1    jpeg_0.1-8.1        
 [22] colorspace_1.4-1     rvest_0.3.5          haven_2.2.0         
 [25] xfun_0.12            callr_3.4.2          crayon_1.3.4        
 [28] zoo_1.8-7            glue_1.3.2           gtable_0.3.0        
 [31] webshot_0.5.2        pkgbuild_1.0.6       rstan_2.19.3        
 [34] abind_1.4-5          scales_1.1.0         mvtnorm_1.1-0       
 [37] DBI_1.1.0            miniUI_0.1.1.1       viridisLite_0.3.0   
 [40] htmlTable_1.13.3     HDInterval_0.2.0     foreign_0.8-76      
 [43] stats4_3.6.2         StanHeaders_2.21.0-1 DT_0.13             
 [46] htmlwidgets_1.5.1    httr_1.4.1           threejs_0.3.3       
 [49] arrayhelpers_1.1-0   RColorBrewer_1.1-2   ellipsis_0.3.0      
 [52] acepack_1.4.1        farver_2.0.3         pkgconfig_2.0.3     
 [55] loo_2.2.0            nnet_7.3-13          dbplyr_1.4.2        
 [58] labeling_0.3         tidyselect_1.0.0     rlang_0.4.5         
 [61] reshape2_1.4.3       later_1.0.0          cellranger_1.1.0    
 [64] munsell_0.5.0        tools_3.6.2          cli_2.0.2           
 [67] generics_0.0.2       ggridges_0.5.2       evaluate_0.14       
 [70] fastmap_1.0.1        yaml_2.2.1           fs_1.3.2            
 [73] processx_3.4.2       nlme_3.1-145         mime_0.9            
 [76] xml2_1.2.5           compiler_3.6.2       bayesplot_1.7.1     
 [79] shinythemes_1.1.2    rstudioapi_0.11      reprex_0.3.0        
 [82] stringi_1.4.6        highr_0.8            ps_1.3.2            
 [85] Brobdingnag_1.2-6    nloptr_1.2.2.1       markdown_1.1        
 [88] shinyjs_1.1          vctrs_0.2.4          pillar_1.4.3        
 [91] lifecycle_0.2.0      bridgesampling_1.0-0 data.table_1.12.8   
 [94] httpuv_1.5.2         R6_2.4.1             latticeExtra_0.6-29 
 [97] bookdown_0.18        promises_1.1.0       boot_1.3-24         
[100] colourpicker_1.0     MASS_7.3-51.5        gtools_3.8.1        
[103] assertthat_0.2.1     withr_2.1.2          shinystan_2.5.0     
[106] parallel_3.6.2       hms_0.5.3            rpart_4.1-15        
[109] coda_0.19-3          minqa_1.2.4          snakecase_0.11.0    
[112] rmarkdown_2.1        shiny_1.4.0.2        lubridate_1.7.4     
[115] base64enc_0.1-3      dygraphs_1.1.1.6